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Abstract. Relativistic modelling of rotational motion of extended bodies represents one of 
the most complicated problems of Applied Relativity. The relativistic reference systems of IAU 
(2000) give a suitable theoretical framework for such a modelling. Recent developments in the 
post-Newtonian theory of Earth rotation in the limit of rigidly rotating multipoles are reported 
below. All components of the theory are summarized and the results are demonstrated. The 
experience with the relativistic Earth rotation theory can be directly applied to model the 
rotational motion of other celestial bodies. The high-precision theories of rotation of the Moon, 
— ■ Mars and Mercury can be expected to be of interest in the near future. 
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1. Earth rotation and relativity 

Earth rotation is the only astronomical phenomenon which is observed with very high 
accuracy, but is traditionally modelled in a Newtonian way. Although a number of at- 
tempts to esti mate and calculate the relativistic effects in Earth rotation have been un- 
■ dertaken (e.g.. iBizouard et al\ (1992); Brumberg k. Simor] ( 2007 ) and reference therein) 



no consistent theory has appeared until now. As a result the calculations of different 
authors substantially differ from each other. Even the way geodetic precession/nutation 
is usually taken into account is just a first-order approximation and is not fully consistent 
with relativity. On the other hand, the relativistic effects in Earth's rotation are rela- 
tively large. For example, the geodetic precession (1.9" per century) is about 3 x 10~ 4 
of general precession. The geodetic nutation (up to 200 /xas) is 200 times larger than the 
goal accuracy of modern theories of Earth rotation. One more reason to carefully inves- 
tigate relativistic effects in Earth rotation is the fact that the geodynamical observations 
yield important tests of general relativity (e.g., the best estimate of the PPN 7 using 
large range of angular distances from the Sun comes from geodetic VLBI data) and it is 
dangerous to risk that these tests are biased because of a relativistically flawed theory of 
y—i • Earth rotation. 

Early attempts to model rotati onal motion of the Earth in a relativistic framework 
(see, for example, Brumbere 1972) made use of only one relativistic references system to 



. . \ describe both rotational and translational motions. That reference system was usually 
chosen to be quite similar to the BCRS. This resulted in a mathematically correct, but 
k>( I physically inadequate coordinate picture of rotational motion. For example, from that 
coordinate picture a prediction of seasonal LOD- variations with an amplitude of about 

C$ I 75 microseconds has been put forward. 

At the end of the 1980s a better reference system for modelling of Earth rotation 
has been constructed, that after a number of modifications and improvements has been 
adopted as GCRS in the IAU 2000 Resolutions. The OCRS implements the Einstein's 
equivalence principle and represents a reference system in which the gravitational influ- 
ence of external matter (the Moon, the Sun, planets, etc.) is reduced to tidal potentials. 
Thus, for physical phenomena occurring in the vicinity of the Earth the GCRS repre- 
sents a reference system, the coordinates of which are, in a sense, as close as possible to 



1 



2 



S.A. Klioner, E. Gerlach, M.H. Soffcl 



measurable quantities. This substantially simplifies the interpretation of the coordinate 
description of physical phenomena localized in the vicinity of the Earth. One important 
application of the GCRS is modelling of Earth rotation. The price to pay when using 
GCRS is that one should deal not only with one relativistic reference system, but with 
several reference systems, the most important of which are the BCRS and the GCRS. 
This makes it necessary to clearly and carefully distinguish between parameters and 
quantities defined in the GCRS and those defined in the BCRS. 



2. Relativistic equations of Earth rotation 

The model which is used in this investigation was discussed and published bvlKlio" ner et al 
( 20011 ). Let us, however, repeat these equations once again not going into physical de- 
tails. The po st-Newtonian equati ons of motion (omitting numerically negligible terms as 
explained in lKlioner et ~d\ (|200lh ) read 



d ( riah 

dT ^ 



Eabc MbL GcL 



L a (c, w ,n iner ), 



(2.1) 



where T = TCG, C = C ab is the post-Newtonian tensor of inertia an d u> — Lu a is 
the angular velocity of the post-Newtonian Tisserand axes (|Klionerlll996h . M L are the 
multipole moments of the Earth's gravitational field defined in the GCRS, Gl are the 
multipole moments of the external tidal gravito-electric field in the GCRS. In the simplest 
situation (a number of mass monopoles) Gl are explicitly given by Eqs. (19)-(23) of 
Klioner et all (|200ll) . 

The additional torque L a depends on C, u>, as well as on the angular velocity Jlinor 
describing the relativistic precessions (geodet ic, Lense-Thirr i ng an d Thomas precessions) . 



The definition of r2i Ilcr can be found, e.g., in lKlioner et all ([20011 ). A detailed discussi on 



of L a , its structure and consequences will be published elsewh ere ( Klioner et al. 20091) . 

The model of rigidly rotating multipoles (jKlioner et aLll200lh represents a set of formal 
mathematical assumptions that make the general mathematical structure of Eqs. (12.11) 
similar to that of the Newtonian equations of rotation of a rigid body: 
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P aibl M blb2 ... bl , M blb2 ... bl = const, 



(2.2) 

I ^ 2, (2.3) 



where the orthogonal matrix P ab (T) is assumed to be related to the angular velocity uj a 
used in (12.11) as 



u a = \e abc P db {T)-^P dc {T). (2.4) 

The meaning of these assumptions is that both the tensor of inertia C ab and the multipole 
moments of the Earth's gravitational field Ml are "rotating rigidly" and that their rigid 
rotation is described by the same angular velocity w a that appears in the post- Newtonian 
equations of rotational motion. It means that in a reference system obtained from the 
GCRS by a time-dependent rotation of spatial axes both the tensor of inertia and the 
multipole moments of the Earth's gravitational field are constant. 

No acceptable definition of a physically rigid body exists in General Relativity. The 
model of rigidly rotating multipoles represent a minimal set of assumptions that allows 
one to develop the post-Newtonian theory of rotation in the same manner as one usually 
does within Newtonian theory for rigid bodies. In the model of rigidly rotating multipoles 
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only those properties of Newtonian rigid bodies are saved which are indeed necessary for 
the theory of rotation. For example, no assumption on local physical properties ( "local 
rigidity") is made. It has not been proved as a theorem, but it is rather probable that 
no physical body can satisfy assumptions (|2.2I) ~ (|2.4I) . The assumptions of the model of 
rigidly rotating multipoles will be relaxed in a later stage of the work when non-rigid 
effects are discussed. 



3. Post-Newtonian equations of rotational motions in numerical 
computations 

Looking at the post-Newtonian equations of motion (12.1I) - (|2.4I) one can formulate 
several problems to be solved before the equations can be used in numerical calculations: 

A. How to parametrize the matrix P ab ? 

B. How to compute Ml from the standard models of the Earth's gravity field? 

C. How to compute Gl from a solar system ephemeris? 

D. How to compute the torque e a bc M^l G c l out of Ml and Gl? 

E. How to deal with different time scales (TCG, TCB, TT, TDB) appearing in the 
equations of motion, solar system ephemerides, used models of Earth gravity, etc.? 

F. How to treat the relativistic scaling of various parameters when using TDB and/or 
TT instead of TCB and TCG? 

G. How to find relativistically meaningful numerical values for the initial conditions 
and various parameters? 

These questions are discussed below. 



4. Relativistic definitions of the angles 

One of the tricky points is the definition of the angles desc r ibing the Ea rth orientation 
in the relativistic framework. Exactly as in IBretagnon et all (|l997l . Il998h wc first define 



the rotated BCRS coordinates (x, y, z) by two constant rotations of the BCRS as realized 
by the JPL's DE403: 

y I = i?.,.(2r2(i'21.4rt02x") /?,|-0.iir,2!)4".) | y ) . (4.1.) 

DE403 

Then the IAU 2000 transformations between BCRS and GCRS are applied to the coordi- 
nates (t,x,y,z), t being TCB, to get the corresponding GCRS coordinates (T, X, Y, Z). 
The spatial coordinates (X, Y, Z) are then rotated by the time-dependent matrix P 4 - 7 to 
get the spatial coordinates of the terrestrial reference system (£, rj, £). The matrix P % i is 
then represented as a product of three orthogonal matrices: 

I J = RM Rx(") R.W (yJ. (4.2) 

The angles 0, ip and to are used to parametrize the orthogonal matrix P ab and therefore, 
to define the orientation of the Earth orientation in the GCRS. The me aning of the 



terrestrial system (£, 77, £) here is the same as in IBretagnon et all (|1997| ): this is the 
reference system in which we define the harmonic expansion of the gravitational field 
with the standard values of potential coefficients Ci m and Si m . 
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5. STF model of the torque 

The relativistic torque requires computations with symmetric and trace-free cartesian 
(STF) tensors Ml and Gl- For this project special numerical algorithms for numerical 
calculations have been developed. The detailed algorithms and their derivation will be 
published elsewhere. Let us give here only the most important formulas. For each / the 
component D a = e a bc Mj,l—i GcL-i of the torque in the right-hand side of Eq. (12.11) can 



be computed as (A t =AIttU/(21 + 1)!!, aj m = + 1) - m(m + 1) ) 

/ 1— i 

Di = t E <n Mtai + 1+^UC) 

1 \m=0 

+ E a tm (ML GL+i ~ Mf; m+1 g[ m )\ (5.1) 

m— 1 / 
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! \m=0 

+ E a ^ (--^hn S{ tm +i + M[ m+1 GL)) , (5.2) 
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^ = X E m ( M L $L - Mt GL) • (5-3) 



The coefficie nts QjL and ffL, cha racterizing the tidal field can be computed from Eqs. 



(19)-(23) of iKlioner et al\ (|2001l ) as explicit functions of the parameters of the solar 
system bodies: their masses, positions, velocities and accelerations. A Fortran code to 
compute GL and Q\ m for I < 7 and ^ m ^ I has been generated automatically with 
a specially written software package for Mathematica. It is possible to develop a sort 
of recursive algorithm to compute GL and G\ m for any / similar to the corresponding 
algorithms for, e.g., Legendre polynomials. 

The coefficients M.f m and M\ m characterizing the gravitational field of the Earth can 
be computed as 



11 ( 4tt 



ML^-ir\j^ {4h^T MERlEC ^ l<m<h (5 - 5) 

M\ m = (-l)" l+1 \ (4^ 7^) V2 M B R l E S lm , l^m^l, (5.6) 



2 (2Z-1)!! \ 2l + l {I -my., 

where Me is the mass of the Earth, Re its radius, C/ m and Si m are usual potential coef- 
ficients of the Earth's gravitational field. If only Newtonian terms are considered in the 
torque this formulation with STF tensors is fully equiva l ent to the classical formulation 
with Legendre polynomials (e.g., Bretagnon et al. 1997 . 1998h . If the relativistic terms 



are taken in account, the only known way to express the torque is that with STF tensors. 



6. Time transformations 

An important aspect of relativistic Earth rotation theory is the treatment of different 
relativistic time scales. The transformation between TDB and TT at the geocenter (all 
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the transformations in this Section are meant t o be "e valuated at the geocenter") are 
computed along the lines of Section 3 of Klioner ( 2008bl ). Namely, 



TT = TDB +ATDB(TDB), 
TDB = TT -ATT(TT), 
TCG = TCB +ATCB(TCB), 
TCB = TCG -ATCG(TCG), 



so that 



dATDB 
dTDB 

^4tdb 

Btdb 
cZATT 



^4tdb + -Bxdb 

Lb — Lq 
1 — Lb 
1-L g 



dATCB 
tfTCB ' 



I -L B 

Lb — Lq 
I -La' 
I -L B 



^4tdb + 1, 
dATCG 



dTCG ' 



1 



= 1 - A 



TT, 



dTT 

dATCB 
dTCB 
dATCG 

dTCG " 1 + F(TCG - ATCG) : 



F(TCB) = \ a(TCB) + \ /3(TCB), 
cr c 4 

F(TCG - ATCG) 



(6.1) 
(6.2) 
(6.3) 
(6.4) 



(6.5) 
(6.6) 
(6.7) 
(6.8) 
(6.9) 
(6.10) 
(6.11) 
(6.12) 



Klioner 



2008H) and Eq. 
Klionerl (l2008bf ). Clearly, 



where the functions a and /? are given by Eqs. (3.3)-(3-4) of 
(|6.12p represents a computational improvement of Eq. (3.8) of 
the derivatives dATCB/dTCB and d ATCG jdTCG must be expressed as functions of 
TDB and TT, respectively, when used in ([63 ]) -([6T8 ]) . 

The differential equations for ATDB and ATT are first integrated numerically for the 
whole range of the used solar system ephemeris (any ephemeris with DE-likc interface can 
be used with the code) . The initial conditions for ATDB and ATT are chosen according to 
the IAU 2006 Resolution defining TDB: for JD TT = 2443144.5003725 one has JL>tdb = 
2443144.5003725-6.55 x 10 5 /86400 and vice versa. The results of the integrations for the 
pairs ATDB and dATCB /dTCB, and ATT and dATT/dTT are stored with a selected 
step in the corresponding time variable (TDB for ATDB and its derivative, and TT 
for ATT and its derivative). A cubic spline on the equidistant grid is then constructed 
for each of these 4 quantities. The accuracy of the spline representation is automatically 
estimated using additional data points computed during the numerical integration. These 
additional data points lie between the grid points used for the spline and are only used 
to control the accuracy of the spline. The splines precomputed and validated in this 
way are stored in files and read in by the main code upon request. These splines are 
directly used for time transformation during the numerical integrations of Earth rotation. 
Although this spline representation requires significantly more stored coefficients than, for 
example, a representation with Chebyshev polynomials with the same accuracy, the spline 
representation has been chosen because of its extremely high computational efficiency. 
More sophisticated representations may be implemented in future versions of the code. 
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7. Relativistic scaling of parameters 

Obviously, there are two classes of quantities entering Eqs. (12.1I) - (|2.4I) that are defined 
in the BCRS and GCRS and, therefore, naturally parametrized by TCB and TCG, 
respectively. It is important to realize that the post-Newtonian equations of motion are 
only valid if non-scaled time scales TCG and TCB are used. If TT and/or TDB are 
needed, the equations should be changed correspondingly. 

The relevant quantities defined in the GCRS and parametrized by TCG are: (1) the 
orthogonal matrix P ab and quantities related to that matrix: angular velocity oj a and 
corresponding Euler angles ip, ip and uj; (2) the tensor of inertia C ab ; (3) the multipole 
moment of Earth's gravitational field Ml- In principle, (a) Gl and (b) Q? neI are also 
defined in the GCRS and parametrized by TCG, but these quantities are computed 
using positions x^_, velocities va and accelerations a a of solar system bodies. The orbital 
motion of solar system bodies are modelled in BCRS and parametrized by TCB or TDB. 
The definition of Gl is conceived in such a way that positions, velocities and accelerations 
of solar system bodies in the BCRS should be taken at the moment of TCB corresponding 
to the required moment of TCG with spatial location taken at the geocenter. Let us recall 
that the transformation between TCB and TCG is a 4-dimensional one and require the 
spatial location of an event to be known. 

In all theoretical works aimed to derive and/or analyze the rotational equations of 
motion in the GCRS one uses TCG as coordinate time scale parametrizing the equa- 
tions. Although the natural time variable for the equations of Earth rotation is TCG, 
in principle, using a corresponding re-parametrization any time scale (including TCG, 
TT, TCB and TDB) can be used as independent time variable. Thus, simple rescaling of 
the first and second derivatives of the angles entering the equations of rotational motion 
should be applied to use TT instead of TCG: 

= (1-£g)-t^, (7-1) 



dTCG 
dH 



(1 



dTT 

„ d 2 9 



(7.2) 



dTCG 2 dTT 2 ' 
where is any of the angles (p, ip and uj used in the equations of motion to parametrize 
the orientation of the Earth. If TDB is used as independent variable the corresponding 
formulas are more complicated: 
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( dTT 



dO 



\ dTDB 



/ dTT 
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\ dTDB 



dTDB 



( dTT 




3 d 2 TT 




d6 


I dTDB 




dTDB 2 




dTDB ' 



(7.4) 



where the derivatives of TT w.r.t. TDB should be evaluated at the geocenter (i.e., for 
x = xe)- These relations must be substituted into the equations of rotation motion to 
replac e the derivatives of the angles tp, ip and uj w.r.t. TCG as appear e.g., in Eqs. (7)— 
(9) of lBretagnon et al. (1998). It is clear that the parametrization with TDB makes the 
equations more complicated. 

The values of the parameters naturally entering the equations of rotational motion 
must be interpreted as unsealed (TCB-compatible or TCG-compatible) values. If scaled 
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(TT-compatiblc or TDB-compatible) values are used, the scaling m ust be explicitly taken 
into account. The relativistic scaling of parameters read (see, e.g.. lKlionerll2008al) : 



GM 



TT 

A 



(1 



L G )GMj GG , GMj GG - 

fTDB 



GMj CB , 



GM^ = (1 - L B ) GM 



TCB 
A i 



-TCG 



X LL ={1-L G )X 
V TT = V TGG , v TDB =v 
A TT = (1 - Lg)- 1 A TGG 



X TDB 
TCB 



.TDB 



(l-L B )x 



(I -L B ) 



TCB 



,TCB 



(7.5) 
(7.6) 
(7.7) 
(7.8) 



where GMa is the mass parameter of a body, x, v, and a are parameters represents 
spatial coordinates (distances), velocities and accelerations in the BCRS, respectively, 
while X, V, and A are similar quantities in the GCRS. 

Now, considering the source of the numerical values of the parameters used in the 
equations of Earth rotation we can see the following. 

a. The position xa, velocities va, accelerations a a and mass parameters GMa of the 
massive solar system bodies are taken from standard JPL ephemerides and are TDB- 
compatible. 

b. The radius of the Earth comes together with the potential coefficients Ci m and 
S im horn a model of the Earth's gravity field (e.g., GEMT3 was used in SMART). These 
values come from SLR and dedicated techniques like GRACE. GCRS with TT-compatible 
quantities is used to process these data. Therefore, the values of the radius of the Earth is 
TT-compatiblc. Obviously, C; m and Si m have the same values when used with any time 
scale. The mass parameter GMe of the Earth coming with the Earth gravity models is 
also TT-compatible. 

c. From the defini tions of A4fL and given above and formulas for Gl given 



by Eqs. (19)-(23) of iKlioner et ali (|200lh . it is easy to see that the TCG-compatible 
torque F a = Y^ili jf £ abc M^l G c l can be computed using TDB-compatible values of 
mass parameters GAf™ B , positions x A T ' B , velocities u™ B and accelerations ct™ B of 
all external bodies, TDB-compatible value of the mass parameter of the Earth GA/| DB 
and the value of the Earth's radius formally rescaled from TT to TDB as -R™ B = 
(1 — Lb) (1 — Lq)^ 1 -R^ T . Denoting the resulting torque by F^ DB , it can be seen that 
the TCG-compatible value is F^, CG = (1 — Lb)^ 1 -Ftdb- 

d. The values of the Earth's moments of inertia Ai, i — 1, 2, 3 can be represented as 
GAi = GMER 2 E ki, where k% is a factor characterizing the distribution of matter inside 
the Earth. Clearly, the factors ki do not depend on the scaling. Therefore, the moments 
of inertia can be scaled as 

AJ T = (l-L G fAj GG . (7.9) 

The last question is how to interpret the values of the moments of inertia Ai — 
( A, B, C) and the initial c onditions for the angles ip, ip and u) and their derivatives given 
in lBretagnon et al. (1998). Obviously, the initial angles at J2000 are independent of the 



scaling. For the other parameters in question it is not possible to clearly claim if the given 
values are TDB-compatible or TT-compatible. Arguments in favor of both interpreta- 
tions can be given. A rigorous solution here is only possible when all calculations leading 
to these quantities are repeated in the framework of General Relativity. In this paper 
we prefer to interpret the SMART values of Ai, ip, ip an d Cj as being TT-compatible. 
Therefore, if TDB is used as independent variable, the values of the derivatives should 
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be changed accordingly. For any of these angles one has 



da I dTT 



dTDB \ dTDB 



Thus, we have all tools to treat correctly the relativistic scaling of all relevant parameters 
of the Earth rotation theory as well as relativistic time scales. 



8. Geodetic precession and nutation 

In the framework of our model geodetic precession and nutation are taken into account 
in a natural way by including the additional torque that depends on f2° nor in the equations 
of rotational motion: 

L a = s abc C bd Lo d n? nor - A (C ab Of ner ) . (8.1) 
The first term of the additional torque reflects the fact that the GCRS of the IAU 



is defined to be kinematically non-rotating (see ISoffel et al. I l2003f) . The second term 



has been usually hidden by the corresponding re-d efinition of the post-Newtonian spin 
(|Damour. Soffel fe Xulll993t iKlioner fc Soffell l2000h . It can be demonstrated that this 
second term must be explicitly taken into account to maintain the consistency between 
dynamical ly and kinematically non-rotating solutions. Further details will be published 
elsewhere ( Klioner et al. 20091 ). Using the additional torque L a in Eq. (|2.1j) is a rigorous 



way to take geodetic precession/nutation into account. 

The standard way to account for geodetic precession/nutation that was used up to 
now by a number of authors can be described as follows: (1) solve the purely Newto- 
nian equations of rotational motion and consider this solution as a relativistic one in a 
dynamically non-rotating version of the GCRS and (2) add the precomputed geodetic 
precession/nutation to it. The second step is fully correct since the geodetic preces- 
sion/nutation is by definition the rotation between the kinematically and dynamically 
non-rotating versions of the GCRS and it can be precomputed since it is fully indepen- 
dent of the Earth rotation. The inconsistency of the first step comes from the fact that 
in the computation of the Newtonian torque the coordinates of the solar system bod- 
ies are taken from an ephemeris constructed in the BCRS. However, the dynamically 
non-rotating version of the GCRS rotates relative to the BCRS with angular velocity 
fi° ner (T). This means that the BCRS coordinates of solar system bodies should be first 
rotated into "dynamically non-rotating coordinates" and only after that rotation those 
coordinates can be used to compute the Newtonian torque. For this reason this procedure 
does not lead to a correct solution in the kinematically non-rotating GCRS (see Fig. [T]). 
We will call such solutions in this paper "SMART-like kinematical solutions" . 

On the other hand, there are two ways to obtain a correct kinematically non-rotating 
solution: (1) use the torque given by Eq. (|8.ip in the equations of motion, (2) compute 
the geodetic precession/nutation matrix, apply the geodetic precession/nutation to the 
solar system ephemeris, integrate (I2.1j) without L a with the obtained rotated ephemeris 
(the correct solution in a dynamically non-rotating version of the GCRS is obtained 
in this step), apply the geodetic precession/nutation matrix to the solution. We have 
implemented both ways in our code and checked explicitly that they give the same 
solution (to within about 0.001 /Ltas over 150 years). It is interesting to note that the 
rotational matrix of geodetic precession/nutation (that is, the matrix defining a rotation 
with the angular velocity rif ncr ) cannot be parametrized by normal Euler angles. We have 
used therefore the quaternion representation for that matrix. 
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Figure 1. Scheme of the two ways to obtain a kinematically non-rotating solution from a purely 
Newtonian one, and an illustration of the relation between the correct kinematically non-rotating 
solution and SMART-like kinematical solutions. "GP" stands for geodetic precession/nutation. 
Each gray block represents a solution. A solid line means: add precomputed geodetic preces- 
sion/nutation into a solution to get a new one. A dashed line means: recompute a solution with 
indicated change in the torque model. 



9. Overview of the numerical code 

A code in Fortran 95 has been written to integrate the post-Newtonian equations of 
rotational motion numerically. The software is carefully coded to avoid numerical insta- 
bilities and excessive round-off errors. Two numerical integrators with dense output - 
ODEX and Adams-Bashforth-Moulton multistep integrator - can be used for numerical 
integrations. These two integrators can be used to crosscheck each other. The integrations 
are automatically performed in two directions - forwards and backwards - that allows 
one to directly estimate the accuracy of the integration. The code is able to use any type 
of arithmetic available with a given current hardware and compiler. For a number of 
operations, which have been identified as precision-critical, one has the possibility to use 



either the library FMLIB ISmithl (|2001l) for arbitrary-precision arithmetic or the pack- 
age DDFU N that uses tw o double-precision numbers to implement quadrupole-precision 
arithmetic ( Bailey 2005). Our current baseline is to use ODEX with 80 bit arithmetic. 



The estimated errors of numerical integrations after 150 years of integration are below 
0.001 //as. 

Several relativistic features have been incorporated into the code: (1) the full post- 
Newtonian torque using the STF tensor machinery, (2) rigorous treatment of geodetic 
precession/nutation as an additional torque in the equations of motion, (3) rigorous 
treatment of time scales (any of the four time scales - TT, TDB, TCB or TCG - evaluated 
at the geocenter can be used as the independent variable of the equations of motion 
(TCG being physically preferable for this role), (4) correct relativistic scaling of constants 
and parameters. All these "sources of relativistic effects" can be switched on and off 
independently of each other. 

In order to test our code and the STF-tensor formulation of the torque we have 
coded also the c l assical New tonian torque with Legendre polynomials as described by 
Bretagnon et al. ( 1997 . 19981) and integrated our equations for 150 years with these two 



torque algorithms. Maximal deviations between these two integrations were 0.0004 /ias 
for (j>, 0.0001 /zas for ip, and 0.0002 /j,as for u. This demonstrates both the equivalence of 
the two formulations and the correctness of our code. 

We have also repeated the Newtonian dynamical solution of SMART9 7 using the New- 
tonian torque, the J PL eph e meris DE403 and the same initial values as in lBretagnon et al 



19981 ). Jean-Louis ISimon (2007) has provided us with the unpublished full version of 



SMART97 (involving about 70000 Poisson terms for each of the three angles). We have 
calculated the differences between that full SMART97 series and our numerical integra- 
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Figure 2. Differences (in fj,as) between the published kinematical SMART97 solution and the 
correct kinematically non-rotating solution (with post-Newtonian torques, relativistic scaling 
and time scaled neglected). 



tion over 150 years. Analysis of the results and a comparison to lBretagnon et ali (|1998I ) 



have demonstrated that our integrations reproduce SMART97 within the full accuracy 
of the latter. 



10. Relativistic vs. Newtonian integrations 

We have performed a series of numerical calculations comparing purely Newtonian 
integration with integration where relativistic effects are taken into account. The same 
initial conditions and parameters that we used to reconstruct the SMART97 solution 
were used for all integrations (see below). The results are illustrated on Figs. [2HH1 The 
difference between the kinematical SMART97 solution and the consistent kinematically- 
non-rotating solution obtained as described in Section [8] is shown in Fig. [2] Fig. [3] shows 
the effects of the post-Newtonian torque. The effects of the relativistic scaling and time 
scales are depicted in Fig. @] Finally, Fig. [5] demonstrates the differences between a 
SMART-like kinematical solution and our full post-Newtonian integration. A detailed 
analysis of these results will be done elsewhere. 

To complete the consistent post-Newtonian theory of Earth rotation the parameters 
(first of all, the moments of inertia of the Earth) should be fitted to be consistent with 
the observed precession rate. This task will be discussed and treated in the near future. 



11. Relativistic effects in rotational bodies of other bodies 

The same numerical code can be applied to model the rotational motion of other 
bodies. Especially, high-accuracy models of rotational motion of the Moon, Mercury and 
Mars are of interest because of the planned space missions to Mercury and Mars, and 
the expected improvements of the accuracy of LLR. Most of the changes in the code are 
trivial and concern the numerical values of the constants. One important improvement 
of the code is necessary for the Moon: the figure-figure interaction with the Earth must 
be taken into account. Using the STF approach to compute the torque this task is not 
difficult. 

The relativistic effects in the rotation of Moon, Mars and Mercury may be significantly 
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Figure 3. The effect (in uas) of the post-Newtonian torque. 
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Figure 4. The effect (in uas) of the relativistic scaling and time scales. 



larger than in the rotation of Earth. In Table [T] the amplitudes of geodetic precession 
and nutation are given for several solar system bod ies. One can see the l arge e ffects for 
Mercury and Mars. Besides an early investigation of Bois &: Vokroulickv ( 1995t) suggests 
that the effects of the relativistic torque for the Moon may attain 1 mas. Our approach 
allows one to investigate the rotational motion of the Moon, Mars and Mercury in a 
rigorous relativistic framework. 
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body geodetic precession geodetic nutation 
[" per century] [fias] 



Earth 

Moon 

Mercury 

Venus 

Mars 



1.92 
1.95 
21.43 
4.32 
0.68 



153 
154 
5080 
85 
567 



Table 1. Magnitude of geodetic precession/nutation for various bodies 
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